cat("数学関数\n")

cat("和--sum\n")


cat("算術平均--mean,colMeans\n")

cat("累積--prod\n")

cat("幾何平均(実装)--calgeomean,subgeomean,geomean\n")

#NAデータなしのベクタの幾何平均を計算する
calgeomean<-function(x){
  if(any(is.na(x)))
    stop("Na data can not be used in calgeomean")
  
  if(any(x<=0))
    stop("Geomean can only caculated for positive data.")
  
  #tprod<-prod(x)
  #n<-length(x)
  #return(exp((log(tprod)/n)))
  return(prod(x)^(1/length(x)))
}

#ベクタの幾何平均を計算する
subgeomean<-function(x,na.rm=TRUE){
  
  if(all(!is.na(x))){
    return(calgeomean(x))
  }else{
    if(na.rm){
      x<-x[!is.na(x)]
      return(calgeomean(x))
    }else{
      stop("There is NA in the data, set na.rm=TRUE ?")
    }
  }
}

#データフレーム、マトリクスおよびベクターの幾何平均を計算する
#データフレームやマトリクスに関して、行ごとの幾何平均を計算する
geomean<-function(x,na.rm=TRUE){
    
  if(is.vector(x))
    return(subgeomean(x,na.rm))
  
  if(is.matrix(x)||is.data.frame(x)){
    if(any(is.na(x)) && (!na.rm)){
      stop("There is NA in the data, please set na.rm=TRUE ?")
    }
    return(apply(x,1,subgeomean))
  }
  
  stop("Data is worng for caculating geomean")
}

cat("調和平均(実装)--subharmean,harmean\n")
subharmean<-function(x,na.rm=TRUE){
  if(any(is.na(x))&&(!na.rm)){
    stop("There is NA in the data, please set na.rm=TRUE ?")
  }
  
  x<-x[!is.na(x)]
  
  if(any(x==0)){
    stop("Harm mean can not be caculated for 0.")
  }
  
  return(1/(sum(1/x)/length(x)))
}

harmean<-function(x,na.rm=TRUE){
  if(any(is.na(x))&&(!na.rm)){
    stop("There is NA in the data, please set na.rm=TRUE ?")
  }

  if(is.vector(x)){
    return(subharmean(x,TRUE))
  }
  
  if(is.matrix(x)||is.data.frame(x)){
    return(apply(x,1,subharmean))
  }
  
  stop("Data is worng for caculating geomean")
}


cat("加重平均--weighted.mean\n")
# weightmean<-function(x,w,na.rm=T){
#   if(length(x)!=length(w)){
#     stop("Size of data and weight are differ")
#   }
#   
#   if(any(is.na(x))&&(!na.rm)){
#     stop("There is NA in the data, please set na.rm=TRUE ?")
#   }
#   
#   index<-!(is.na(x)|is.na(w))
#   if(all(index==FALSE)){
#     stop("No meanful data can be caculated.")
#   }
#   
#   x<-x[index]
#   w<-w[index]
#   
#   if(sum(w)==0){
#     stop("Sum of all meanful weight is 0 ")
#   }
#   
#   return(sum(x*w)/sum(w))  
# }

cat("ベクタシフト(実装)--shift_ver\n")
shift_vec<-function(x,k=1,direction="right",loop=FALSE){
  if(!is.vector(x)){
    stop("Only vector can be shifted")
  }
  
  if(is.null(x))
    return(x)
  
  if(!loop){
    if(k>=length(x))
      return(rep(NA,length(x)))
  }
  
  k<-k%%length(x)
  
  if(k==0)
    return(x)
  
  if(direction=="left"){
    if(loop){
      return(c(x[(k+1):length(x)],x[1:k]))
    }else{
      return(c(x[(k+1):length(x)],rep(NA,k)))
    }
  }else{
    if(loop){
      return(c(x[(length(x)-k+1):length(x)],x[1:(length(x)-k)]))
    }else{
      return(c(rep(NA,k),x[1:(length(x)-k)]))
    }
  }
}

cat("誤差関数(実装)--suberf,erf\n")
suberf<-function(x){
  f<-function(t){
    exp(-t^2)
  }
  
  return(2/sqrt(pi)*integrate(f,0,x)$value)
}

erf<-function(x){
  return(sapply(x,suberf))
}

ploterf<-function(){
  x<-seq(-3,3,by=0.01)
  plot(x,erf(x),type="l",col="red",main="誤差関数Erf(x)")
  abline(h=c(-1,-0.5,0,0.5,1),v=c(-3,-2,-1,0,1,2,3),col="grey")
  abline(h=0,v=0,col="black")
}

cat("分散--var,ovar(実装)\n")
ovar<-function(x){
  meanx<-mean(x);
  diff<-sum((x-meanx)^2)
  return(diff/length(x))
}
